Fault-Tolerant Renormalization Group Decoder for Abelian Topological Codes 
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We present a tiiree-dimensional generalization of a renormalization group decoding algorithm for 
topological codes with Abelian anyonic excitations that we introduced for two dimensions in [TJ[2]. 
This 3D implementation extends our previous 2D algorithm by incorporating a failure probability 
of the syndrome measurements, i.e., it enables fault-tolerant decoding. We report a fault-tolerant 
storage threshold of ~ 1.9(4)% for Kitaev's toric code subject to a 3D bit-flip channel (i.e. including 
imperfect syndrome measurements). This number is to be compared with the 2.9% value obtained 
via perfect matching [J- The 3D generalization inherits many properties of the 2D algorithm, 
including a complexity linear in the space-time volume of the memory, which can be parallelized to 
logarithmic time. 



I. INTRODUCTION 

Topological quantum error-correcting codes currently 
stand as some of the most promising implementations of 
quantum memories and computers. Crudely, topologi- 
cal codes are standard quantum error-correcting codes 
with additional geometric constraints: their check oper- 
ators involve only neighbouring spins on a two dimen- 
sional (2D) lattice. As a consequence, they can exhibit 
high fault-tolerant threshold [IH^ with relatively low 
overhead. Some topological codes also support transver- 
sal implementation of Clifford gates [7], which simplifies 
fault-tolerant quantum computation. Lastly, topological 
codes can be efhciently decoded [3 El [S], which is the 
topic of this paper. 

Decoding a quantum code consists in inferring the op- 
timal recovery given a statistical description of the noise 
and an error syndrome — i.e., the measurement outcome 
of check operators which reveal incomplete information 
about the particular error that has affected the system. 
Thus, decoding is a classical statistical inference problem 
involving a very large number of correlated random vari- 
ables. Extremely fast decoding algorithms are required to 
prevent errors from building up in between error correc- 
tion cycles, although some lag-time can be tolerated, e.g., 
by extending ideas from [9] . In [U [2] , we introduced a de- 
coding algorithm for Kitaev's topological code TO] that 
uses renormalization group (RG) techniques from statis- 
tical physics. It's complexity is linear with the number 
of qubits, as compared to the cubic complexity of pre- 
viously known algorithms |11| . Most importantly, it can 
be parallelized to logarithmic time. 

The present paper is a continuation of our work initi- 
ated in [2 [2], and serves many purposes. 1) Our previ- 
ous work focused on error correction in the presence of 
perfect syndrome measurements. When measurements 
are faulty, fault-tolerant techniques are required which 
change the nature of the decoding problem. As we ex- 
plain below, for topological codes, this can be effectively 
described by increasing the lattice dimension by one di- 
mension representing time |8]. Thus, we adapt our RG 
algorithm, initially devised for a 2D lattice, to a 3D fault- 



tolerant setting. -'^ 2) Our algorithm was devised specifi- 
cally for Kitaev's topological code. Because all 2D stabi- 
lizer codes are locally equivalent to multiple copies of Ki- 
taev's code [13], our RG algorithm can be used with any 
such code. However, this requires determining the local 
mapping that realizes this equivalence, and transforming 
the local noise model accordingly, which can in principle 
affect the decoder's performances. Here, we describe our 
methods in physical terms that are directly applicable to 
any code that supports Abelian anyons [101 [T5HTB] , not 
restricted to stabilizer codes. We have implemented a 
special case of this generalization in [16' for the Z^ quan- 
tum double model. 3) Our previous publications on this 
topic focused on applications, giving only a high level 
description of the actual algorithm. Here, we provide a 
complete detailed description of the structure of the al- 
gorithm, which should be sufficient for anyone interested 
in implementing it. 



The rest of the paper is organized as follows. In the 
next section, we provide a heuristic physical description 
of the algorithm in terms of localized Abelian anyons. 
This section should provide a good physical intuition of 
the different components of the algorithm. This is first 
done assuming perfect syndrome measurements, and in 
the last subsection we explain how the problem is mod- 
ified in the presence of faulty errors, following jSj. Sec- 
tion [Hllrcvisits all the concepts introduced heuristically 
in Sec.|II|for the special case of Kitaev's topological code, 
using an algebraic formalism closely related to the actual 
implementation of the algorithm. Section |IV| presents 
our numerical experiments, and we conclude in Sec. W\ 
with possible extensions and relations to other methods. 
Appendix [X] det ails our mathematical notation for prob- 
ability distributions over the n-qubit Pauli group. 



^ Note that we have used our algorithm in a fault-tolerant setting 
in 1121 . but did not provide any details of the implementation. 
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FIG. 1: (a) A 2D topological code is cut into unit cells a, /?, 
... Gauge lines representing the non-trivial cycles (solid red 
lines) are chosen arbitrarily. Computing the flow of charge 
through the gauge lines is equivalent to decoding, (b) Each 
region has four boundaries that we label north (A^), east (_E), 
south (S) and west iW). 



II. HEURISTIC PHYSICAL DESCRIPTION 
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FIG. 2; Structure of the RG cells. A unit cell is composed 
of four regions (unit cells of the previous RG iteration). In 
each unit cell (red square), the charge of only three of the 
four regions is measured (green squares); the south-east cor- 
ner is not measured, leaving the total charge of the unit cell 
undetermined. This missing measurement is replaced at the 
following RG iteration by a measurement of the entire unit 
cell (red square) , which is now a region of a renormalized unit 
cell (blue square). Note that this modification of the charge 
measurement does not need to be implemented physically, it 
only reflects a change in bookkeeping. 



In this Section, we provide a heuristic physical descrip- 
tion of the problem of interest, and of the numerical tools 
we have developed to solve it. A more detailed mathe- 
matical description is presented in Sec. |III[ 



decoding problem consists in determining the net flow of 
charge, or current, across these two gauge lines. 



A. Decoding problem 

Consider a 2D sheet of topological matter supporting 
Abelian anyons. For simplicity, suppose that the system 
has periodic boundary conditions, so it forms a torus. 
The information is encoded in the degenerate ground 
state of the system. Excitations above the ground state 
manifold are localized Abelian anyons — they carry con- 
served charges {a, b,c, . . .} that obey "deterministic" fu- 
sion rules, e.g. a x b = c. The information in the ground 
state can be modified by creating a particle-antiparticle 
pair {a, a), dragging one of the particle around a topo- 
logically non-trivial cycle, and fusing it with its original 
partner a x a = 1. 

In the presence of errors, such a process could occur 
spontaneously. For instance, the creation of a particle- 
antiparticle pair could result from a thermal fluctuation. 
Once created, additional errors could cause the parti- 
cles to diffuse on the sheet. To prevent corruption of the 
memory, we must therefore keep track of the homology of 
the particles' world-lines. Periodic measurements of the 
particles' location yield partial information about their 
trajectories, and the decoding problem becomes one of 
statistical inference: it sets to determine the most likely 
homology of the particles' world-lines given two consec- 
utive snapshots of their locations. Concretely, we can 
arbitrarily choose two gauge lines representing the two 
non-trivial cycles of the torus [c.f. Fig. l(a)|, and the 



B. RG algorithm 

In [U H], we proposed a renormalization group tech- 
nique to tackle this problem. First, we break the lattice 
into 2x2 sublattices, or "unit cells", as illustrated on 
Fig. 1(a) Given a microscopic noise model, we can com- 



pute the probability for the value of the current across 
each of the four walls [North, South, East, West, c.f. 
Fig. l(b)| of each cell, conditioned on the charge config- 
uration observed inside this cell. This produces a prob- 
ability distribution VaiNa, Ea, Sa,Wa) for each cell a, 
where Na , E^ , Sa , Wa take values representing the pos- 
sible currents.^ 

Concretely, the presence of a charge, say, in the north- 
east corner of the unit cell would lead to the assignment 
of a probability 0{p) to a current through the northern or 
eastern walls, and a probability C(p^) for the southern or 



^ To specify the mathematical structure of the current variables, 
we can choose a minimal set {ai, 02, . . . flfc} of A: "elementary" 
charges that generate all other charges under fusion. Then, any 
charge can be written as a"^ X Oj ^ X . . . a^*" , or more succinctly 
represented by the vector (01, a2, ■ • ■ , Ofe) S Z**! X Z''^ x . . . Z''* 
where hj is the order of charge aj, meaning that hj copies 
of aj always fuse to the identity. Then, the current variables 
Na, Sa, Ea, and Wa each take value in Z'^i X Z''2 x . . . Z'^fc . 
In the case of the toric code for instance, there are two elementary 
charges, e and m, and their order is 2 since exe = mxm = l. 



western walls, reflecting the fact that the first two cases 
require only one error process while the second two cases 
require two error processes. Here, p represents the prob- 
ability of an error process such as particle creation, anni- 
hilation, or displacement. The big-O hides multiplicative 
factors accounting for the distinct error processes result- 
ing in the same currents, as well as higher order processes. 
In any case, these probabilities can be computed exactly 
given an underlying local noise model. 

After having computed these current probability dis- 
tributions for every cell, we merge groups of four neigh- 
bouring unit cells into renormalized cells (c.f. Fig.^ and 
iterate the procedure: we sum over all the bulk processes 
that lead to a given current across each of the four renor- 
malized boundaries of each cell. This is done as explained 
above, except that the error probability p is not uniform 
on the lattice, but is given by the current variables of 
the previous RG iteration. By successive iteration, (and 
assuming for simplicity that the lattice linear dimension 
is a power of 2) we arrive at a situation where the North- 
ern and Western walls actually correspond to the gauge 
line representing the non-trivial cycles of the torus. De- 
termining the current across these walls is equivalent to 
decoding, as explained above. 

The difficulty with the procedure we outlined above 
is that charge conservation imposes strong correlations 
between the current variables, so their exact joint prob- 
ability cannot be computed efficiently. To see this, note 
that the current variables are subject to two constraints, 
(a) The sum of the current entering a cell must be equal 
to the total charge inside the region. This leads to a con- 
servation equation N^ + S^ + E^ + W^ — Ca for each 
cell a, where Cq. is the total charge contained in a, and is 
known from observation (error syndrome), (b) The cur- 
rents associated to juxtaposed walls of neighbouring cells 
must be equal and opposite, e.g. S^ = —Ng when 5 is 
the cell directly to the south of cell a, see Fig. 1(a) This 
simply follows from the fact that, e.g. Sa and Ng are 
actually associated to the same physical boundary. Con- 
straints (a) correlate the variables of a given cell while 
constraint (b) correlate variables between different cells, 
so the distribution is globally correlated. 



Thus, approximations are required to solve this prob- 
lem efficiently, as we now explain. First, just as a matter 
of bookkeeping, each cell stores only the random vari- 
ables associated to its northern and western walls, the 
other ones are redundant from constraint (b). This does 
not affect the correlated nature of the problem however 
since (a) becomes Na-\-Wa—Ns—Wp 
(b) now says that e.g. Va{N, 



Ca [c.f. Fig. 1(a) 

a,Wa,Ns,Wp) 



and (b) now says that e.g. 'Pa{Na-,Wa,Ns,Wi3) and 
Vp{Np, Wi3,N^, W~f) must be the marginals of one global 
distribution V{Na,Wa.,Ns,Np,Wi3,N^,W^). To sim- 
plify the problem, we relax this condition to a "mean- 
field" condition, demanding that the two distributions 
yield the same marginals along the wall they share, i.e. 
VaiWp) = 'Pi3{Wp), where the marginals are defined the 



usual way 

Va{Wp)^ Y, Va{Na,Wa,Ns,W^) (1) 

'Pp{W0)= Y. MNp,Wp,N,,W,). (2) 

Ni3,N,,W., 

These mean-field conditions are enforced heuristically us- 
ing belief propagation ^7] . 

Since mean-field approximations are not reliable in 
strongly correlated systems, we make one more modi- 
fication to the problem. Charge conservation imposes a 
hard constraint (a) to the current variables, which is un- 
likely to ever be fulfilled in a mean-field approximation. 
To circumvent this problem, we let the charge ca inside 
each cell fluctuate, i.e., we treat it as a random variable. 
To describe this procedure, recall that each unit cell is 
composed of a collection of four regions (i.e. unit cells 
of the previous RG iteration). Measuring the charge dis- 
tribution inside the unit cell amounts to measuring the 
total charge in each of these regions, which clearly fixes 
the total charge of the unit cell. In the modified proce- 
dure, we measure the charge of all but one of the regions, 
say the south-east region. As a consequence, the total 
charge of the unit cell is undetermined, which relaxes the 
constraints on the current variables as desired. This pro- 
cedure is illustrated on Fig. [2] The charge of the unit cell 
is only fixed at the following RG iteration. 



C. Fault-tolerant decoding 

Our description of the problem so far assumes that 
the charge measurements are perfect. A realistic noise 
model would also include faulty measurements, i.e. every 
charge measurement has some probability of reporting 
the wrong charge. To alleviate this problem, measure- 
ments can be repeated in time. A different outcome be- 
tween two consecutive measurements can then be caused 
either from an actual error having occurred in the time 
between the measurements — e.g. a particle has moved in 
this region — or by an error in one of the two measure- 
ments. 

Consider the space-time cube enclosed between two 
consecutive local charge measurements (c.f. Fig. jsl. We 
can associate a topological charge to this cube equal to 
the difference between the charges revealed by the two 
measurements enclosing it. If the charge of a cube is 
non-trivial, it means that the two consecutive measure- 
ments did not yield the same result. As explained above, 
this could be caused by a "space-like error" taking place 
between the two measurements, or a "time-like error" af- 
fecting the measurements themselves, see Fig. [3j In any 
case, the total current across the six walls of the cube 
must be equal to the charge of the cube. We then see [8] 
that the decoding problem becomes that of determining 
the world-line homology of the particles in space-time. 
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FIG. 3: Space-time diagram of the fault-tolerant error- 
correction procedure; time flows vertically (a) A space-like 
error is an error that affects a qubit in between two mea- 
surements. It creates excitations in the two cubic cells with 
which it overlaps, (b) A time-like error is caused by a faulty 
measurement. It creates excitations in the two cubic cells 
separated by that measurement. 



Thus, the fault-tolerant decoding problem differs from 
the decoding problem with perfect measurements only in 
respect of the lattice dimension. Hence, the RG decoding 
algorithm outUned above can be applied directly. 



III. FORMAL DESCRIPTION FOR KITAEV'S 
TORIC CODE 
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FIG. 4: a) One site, Aij, of the square lattice. A, on 
which is deflned KTC. Qubits, Aij^ and Ai,j,i, live on the 
edges and are associated to sites with the convention de- 
picted, b) Site operator j4i J — Xij^HXij^vXij-i^nXi-ij^v- 
Blue strings represent X operators, c) Plaquette operator 
Bij — Zi^j^nZij+i^vZi+ij^HZij^v- Green strings represent 
Z operators. 



two global constraints Yi- ^i.j 



11 and a, i?. 



U. 



This implies that two logical qubits are encoded in the 
codespace. 

The logical X and Z operators acting on the encoded 
qubits are non-trivial homological cycles (i.e loops around 
the torus) of X operators on the dual lattice and Z op- 
erators on the direct lattice. We arbitrarily choose the 
bare logical operators to be 

Zq = _[J_^oj\H Zi = J^|^i,o,y (3) 

j i 



In this Section, we describe more rigorously the con- 
cepts introduced in the previous Section for the special 
case of Kitaev's toric code (KTC). We begin with the 
2D scenario as it is technically simpler, yet conceptually 
equivalent to 3D. The system is a f x £ square lattice. A, 
with periodic boundary conditions. We assume that i is 
a integer power of 2. Each site A^ j (0 < i,j < i) holds 
two qubits, Aij^a (ck € {H, V}, where H and V stand for 
horizontal and vertical, respectively). The KTC on the 
torus is a stabilizer code [18] and we assume familiarity 
with this class of codes. 



A. Model 

The stabilizer group of KTC is generated by two 
types of operators. On every site, Aj ,, define a site 



i,j,HXijyXij-i^H 



X. 



ij,v, 



and on 



operator, jj^j 



B,- 1 — 



(see Fig. |4|. Let Sg 



operator, Aij — A"i 

every plaquette, define a plaquette 

Zi,j,H Zij + iy Zi^ljH ZijY 

{Aij,Bij} be the set of all plaquette and site op- 
erators. Note that it is invariant under translation. 
The codespace is defined to be the simultaneous -1-1 
eigenspace of all the stabilizer operators. Equivalcntly, 
we can define the Hamiltonian H ~ — '^oes Q^ ^'^'^ ^^^ 
codespace is the degenerate ground space of H. There 
are n = 2(^ qubits on the lattice but only 2i^ — 2 inde- 
pendent generators, i.e. Sg is overcomplete. Indeed, one 
can easily verify that the stabilizer generators obey the 



These correspond to the gauge lines introduced in the 
previous section, c.f. Fig. |l(a) 



Errors are modeled by random Pauli operators affect- 
ing the qubits. A Pauli operator will in general anti- 
commute with a subset of the elements of Sg, causing 
their eigenvalues to flip from +1 in the codespace to - 
1. An element of Sg with -1 eigenvalue corresponds to a 
local excitation, an Abelian anyon. We refer to a plaque- 
tte excitation as a magnetic flux and to a site excitation 
as an electric charge. It is useful to associate binary 
matrices, a^j- (b^j) to an excitation configuration, with 
entries if the eigenvalue of Aij (Bij) is +1 and entries 
1 otherwise. Thus, the excitation configuration associ- 
ated to the product of two errors is the binary sum of 
their respective excitation configurations — the two dis- 
tinct topological charges are their own inverse. 

Since the Pauli operator Xi j u anti-commutes with 
plaquettes {i — l,j) and (i,j), we see that X operators 
can create a pair of magnetic fiuxes, move a magnetic 
fiux to a neighbouring plaquette, and annihilate a pair 
of neighbouring magnetic fluxes. The Z Pauli operator 
plays an equivalent role for electric charges. Thus, the 
microscopic noise model describing the dynamics of the 
anyons can be specified by a memoryless Pauli channel 
'P^.j,a[Q), Q e {I,X,Z,Y = iXZ}—i.e., a probability 
distribution over the four Pauli operators for each qubit 
of the lattice. In this model, the errors E affecting the 
system are thus elements of the n-qubit Pauli group G^- 
The probability of an error E = (^^ ^ Qi,j.a is simply 
given by r{E) = n.,,.„^..j,o(Q»j,a)-' ' 



B. Decoding problem 

When an error E E Q"^ affects the system initially in 
codespace, the task of error-correction is to bring the sys- 
tem back in the codespace by matching every excitations 
in pairs — thus annihilating them all — without changing 
the encoded information. This is realized by applying 
a correction operator, C E G^ ■ If the total operator 
EC is homologically non-trivial, a logical operation will 
be implemented as the system is brought back to the 
codespace, so the information will be corrupted. To be 
successful, the correction C must therefore be homologi- 
cally equivalent to the error E. 

The decoding problem can be formulated in terms of 
this equivalence. Given an error syndrome — i.e., an exci- 
tation configuration — the decoder must find a Pauli op- 
erator that is homologically equivalent to the error that 
has created this syndrome. This is a statistical inference 
problem. One approach to this problem is to find, among 
all errors that are consistent with the observed excita- 
tion configuration, the one with the highest probability. 
When the noise model is independent and uniform, this 
error is simply the lowest weight operator consistent with 
the excitation configuration, where the weight of C is the 
number of non-trivial singic-qubit Pauli operators in C. 
The Perfect Matching Algorithm (PMA) performs this 
task with a 0(^6) complexity [31[H]. 

This turns out not to be the optimal solution however. 
To understand this, let t denote an operator with the 
correct excitation configuration. We suppose that t is 
chosen in some canonical way, so it is in one-to-one cor- 
respondence with excitation configurations. The proba- 
bility that the error E is homologically equivalent to t is 
simply proportional to the sum of the error probability 
'P{Q) over all errors Q equivalent to t. Since the equiv- 
alence relation is generated by elements of the stabilizer 
group S, this is X]se5 ^(^s). On the other hand, t could 
differ from the actual error Ehj a. combination of logical 
operators Eq. (|3|, i.e. a non-trivial cycle. Thus, we can 
use the group generated by the logical operators Eq. ([3| 
to label the equivalence classes of errors. Generalizing 
the above reasoning, the probability that the error E 
is homologically equivalent to tl defines the probability 
associated to the class I G {Xi,Zi) conditioned on the 
excitation configuration (or equivalently conditioned on 
t): 



^ ' ses 



(4) 



where the normalization factor is 7'(t) — J2i s'^i^^^)- 
The optimal decoding consists in choosing the I that max- 
imizes Eq. Q (so the normalization V{t) is not relevant). 
The product tls is a specific Pauli operator and V{tls) 
is the probability of this operator as given by the noise 
model. This computation is intractable because \S\ scales 
exponentially with the system size. 

The type of mathematical manipulation leading to 
Eq. Q will be used extensively by the algorithm and in 




FIG. 5; Left: Choice of a 2 x 2 unit cell used to perfom the RG 
on the KTC. Green disks represent plaquette operators, blue 
squares represent site operators, and edges represent qubits. 
The two generators which are left out are represented by an 
empty square and a circle. Right: The RG yields a renormal- 
ized lattice. The eigenvalue of the renormalized generators 
corresponds to the total charge of the region and the renor- 
malized noise model corresponds to the net flow of charges 
throught the boundaries (eq. [6]). 



the following discussion, so Appendix |A] provides some 
formal background and examples that should be con- 
sulted before reading the next sections. 



C. RG decoding algorithm 

The RG algorithm decomposes the lattice into unit 
cells. We choose them to be 2 x 2 squares enclosing four 
plaquette and four site generators, see Fig. [5J As ex- 
plained in the previous section, the RG decoder requires 
knowledge of all but one of the magnetic and one of the 
electric operators it encloses. By symmetry, we choose 
to leave out the south-east plaquette operator and the 
north-west site operator. As a consequence, the scheme 
will follow our description of Sec. |ll] for the magnetic 
fiuxes, but for the electric charges the lattice is rotated 
by 180° relative to our description of Sec.|lT]. We include 
in the cell all the qubits that participate in the excita- 
tions measured operators, so a cell contains 12 qubits in 
total. Some of the qubits are shared between neighbour- 
ing cells, and this will be responsible for the constraint 
(b) that correlates their current variables. 

To set up calculations, we define the following basis for 
the 12 qubits of the unit cell, see Fig. [6] for qubit labeling 

So — XdX-zXzX^ To = Zo Eo = XeXio Xq = X2X6 

Si — X1X4X5X9 Tl = Zi -El = X7X11 Xi = X5X7 

52 ~ X-iXiXoXj T2 = Z0Z3 E-z = ZoZg, Zo = Z0Z2 

53 — Z0Z1Z3Z4 Tg = X4X7 E3 = ZiZc, Z\ = Z\Z^ 
Sa ~ ZzZzZqZu) T4 = X(, E4 = Xs 

Ss ~ ZaZsZtZii T5 = X7 E^y = Xg (5) 

E(j = Zio 
Er = Zi\ 

The physical interpretations of these operators are the 
following. The stabilizer generators Sj are the six ex- 
citation measurement operators used in the unit cell; 
they are plaquette and site operators. The Tj are the 
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FIG. 6: Two neighbouring unit cells labeled L and R. Each 
shows the labeling of the qubits used in Eq. ^ . Note that 
since these two cells are neighbours, they share qubits. In 
particular, qubits 6 and 10 in cell L are the same as qubits 9 
and 1, respectively, in cell R. 



associated canonical pure errors in the sense that t = 

citation configuration a^ j- b^ ^ inside the cell, without in- 
ducing any magnetic flow through the northern or west- 
ern wall or any electric flow through the southern or east- 
ern wall. The logical operators Xi and Zi monitor re- 
spectively the magnetic current across the north (z — 0) 
and west (i — 1) wall and the electric current across the 
east {i — 0) and south {i — 1) wall. Thus, they corre- 
spond to the current variables used in Sec. In] Lastly, the 
Ej are errors that change the charge of the site and pla- 
quette operators that have been left out of the cell. For 
instance, Eq brings a magnetic flux through the eastern 
wall into the south-east corner. 

An RG iteration takes an excitation configuration and 
a probability distribution over the Pauli group of the 
qubits contained inside the unit cell, and outputs a cur- 
rent probability distribution obtained by summing over 
all equivalent processes that are consistent with the ob- 
served excitation configuration. For example, the oper- 
ator Xq (see Fig. [6] for labeling) is equivalent to the op- 
erator X2X3 as it corresponds to a flow of one magnetic 
flux through the north boundary and into the north- 
west plaquette. This is more directly seen when de- 
composed in the basis Eq. ([5]): Xq = T3XoS'oS'2-B4 and 
X2X3 = T3X0S2 since both decompose into the logi- 
cal operator Xq, which is associated to the magnetic 
current through the northern wall, and the pure error 
T3 which is conjugate to 6*3, the north-west plaque- 
tte. Thus, if a magnetic flux was indeed observed in 
the north-west corner, both of these errors should con- 
tribute to the probability of a magnetic flow through the 
northen boundary. More generally, the probability of a 
current I € {Xi,Zi) conditioned on a charge configura- 

is given 



tion f = 7^"'.3+i7^"»+i'3 7^"»+i.j+i7^''''3 7^''''j+i7^'''+i.j 

U X ^ \j '-± o 



by 



D. Belief propagation 

In the unit cell of Fig. [6] we see that there are eight 
qubits that belong to two unit cells; they are labeled 
0, 1, 6, 7, 8, 9, 10 and 11. For instance, qubit 1 of cell R is 
the same as qubit 10 of cell L immediately to its left. As 
for any other qubits, knowledge of the excitation config- 
uration affects the error probability of these qubits. For 
instance, suppose that the system contains only two mag- 
netic fluxes, one in the north-east corner of cell L and one 
in the north-west corner of cell R. In cell L, the presence 
of the magnetic flux should yield a high probability of X 
error on qubits 2 and 10. In cell R, the presence of the 
magnetic flux should yield a high probability of X error 
on qubits 1 and 0. But since qubits 10 of cell L and 1 
of cell R are actually the same, this charge conflguration 
should globally result in a very peaked probability of an 
X error on that qubit: it sits in between the two magnetic 
fluxes. But locally, given only knowledge of the charge 
configuration on a unique cell, this conclusion cannot be 
reached. 

More generally, given a probability distribution over 
the Pauli group of the unit cell V(tles), we can compute 
the marginal error probability Vq{tles\q) for each qubit q, 
obtained by taking a marginal of Titles) (c.f. App.[A|).^ 
When a qubit is shared between two cells, e.g. such as in 
the above example, its marginal conditional distributions 
obtained from different cells will typically differ. This is 
a manifestation of a violation of constraint (b) described 
in Secjll] As explained there, the exact solution would be 
to demand that the conditional probability distribution 
assigned by each cell be consistent with one global prob- 
ability distribution. Because of global correlations this 
problem is intractable, so we settle for a relaxed condi- 
tion that the marginal probability distributions all agree. 

This condition is enforced by a belief propagation al- 
gorithm. This is a local message passing algorithm 
where messages are exchanged between neighbouring 
cells, there is one message per shared qubit. Initially, 
the outgoing messages at a cell m°"*(p) are equal to 
Vq{tles\q) computed in that cell. These outgoing mes- 
sages are then exchanged between neighbouring cells, and 
become incoming messages, e.g. a cell L sends to its right 
neighbour R the message m^"* that becomes m™^ in i?, 
and receives from R the message m^fj' that becomes to'" 
in L. Subsequent rounds of messages can be calculated 
using the received messages, following the prescription 



^(p)^^<5(t/es|„p) 



V{tles) 
Vq{tles\q) 



n S'(*^es|,0, 



i3'#9 



(7) 



V{l\t) o^^V{tle 



(6) 



where s € (5,) relates topologically equivalent trajecto- 
ries and e € (Ei) changes the value of the undetermined 
charge, and we left out the normalization factor 'P{t). 



■^ The base error prior is independent on each qubit, in which case 
this marginal consists in the noise model on qubit q. But because 
the RG can create a correlated noise model inside a unit cell, we 
need this more sophisticated notion of marginal, see App. lAl 



Here, q,q' G {0,1,6,7,8,9,10,11}, tles\q is the restric- 
tion to qubit q of the Pauli operator tles^ and Vq is the 
marginal on qubit q of the noise model as above (c.f. 
App. [a]). BP can be iterated a few times (e.g. three 
rounds) before executing a RG step. The messages are 
used to update the prior error probability, effectively re- 
placing Eq. (|6| by 

Vm^ E E V{tles)\{w}^{tles\q). (8) 

ee(£;o.Bi>se(SoSiS2> q 



E. Fault-tolerant decoding 

The prescription given for the 2D decoding problem 
can be applied relatively straightforwardly to the 3D 
problem arising from fault-tolerant decoding in the pres- 
ence of faulty syndromes. To simplify the description, we 
will assume that there are only bit-flip errors [X errors), 
so we only need to consider magnetic fluxes. The exact 
same method applies to Z errors and electric charges, 
and moreover both type of errors can be considered si- 
multaneously (including Y errors). 

We label by < fc < r the time at which the charge 
measurements are performed, where r is the total dura- 
tion of the computation (e.g. here we typically set t — I 
to obtain a space-time cube). Errors affect the qubits 
in between measurements, and we use the label k for 
an event that occurs in between measurement fc — 1 and 
k. There are now two types of errors to be considered. 
Space-hke errors jy*^ iVij.a ^ ^2) result in the apphcation 

of the Pauli operator E^ = Y\^ ^ X^^-j-" to the qubits 

between measurements fc — 1 and fc. Time- like errors /x*^ 
(/if^ G Z2) result in inverting the measurement outcome 
at space-time coordinate {i,j,k) when /j,*^ — 1. 

The excitation configuration measured at time fc re- 
sults from the accumulation of space-like errors at times 
prior or equal to fc, plus the measurement errors at 
that time, i.e. b''" — ^^ + conf(nj.,<j, i?'^') — ^j}^ + 

X]fc'<fc conf (ii^''' ). Thus, the difference between two con- 
secutive rounds of measurements is Ab*^ = b'^"-'^ 4- b*^ = 
^fe-i + ^fe + conf(£;'=). In other words, M^^ = a*^J^ + 

^4,] + riljM + "ihy + ^»>ij\H + '^fj-hi.y This defines a 
local space-time cubic check operator. 



In this 3D picture, a Afo; 



*j 



1 plays the role of a 
magnetic flux. Note that each single error — either spatial 
or temporal — creates a pair of fluxes. In particular, the 
set of all errors can be viewed as a product of strings with 
magnetic fluxes located at their endpoints. 

To formalize this description, define a 3D cubic lattice 
of bits, A, with sites ^i,j,k, holding three bits, Aij^k,a 
(a G {H, V, T}) with the convention that bits live on 
faces (see Fig.l?]). The error history, E, on the 3D lattice 
is the combination of all space-like errors rj and time- 
like errors fi, i.e. E,, 
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FIG. 7: Convention chosen for axis and unit cell of tlie 3D 
cubic lattice A. Bits are located on faces. 



to E is A6*^ -. In the following, we consider periodic 
boundary conditions in the spatial dimension to simplify 
the presentation. Then, as in the 2D case, two error 
histories are equivalent if they have the same excitation 
configuration and their product is homologically trivial 
on the three-torus. 

The decoding problem thus stays qualitatively the 
same: find the most likely equivalence class of error his- 
tories consistent with the error syndrome. One subtle 
difference has to do with homologically non-trivial time- 
like loops, which do not carry the same physical meaning 
as space-like homologically non-trivial loops (logical op- 
erations). This difference is only caused by the unphysi- 
cal boundary conditions that were chosen to simplify the 
presentation and the numerical simulations, and would 
not occur with open boundaries. In any case, a time-like 
logical error should not be regarded as a true memory 
corruption. 

As in 2D, perfect matching [3j can be used to solve 
an approximate version of this problem, that consists of 
finding, among all error histories consistent with the ex- 
citation configuration, the one with highest probability. 

The optimal solution however consists in finding the 
most likely equivalent class of errors, and this problem 
can be approximated with RG techniques. The RG de- 
coding has the same logical structure as in 2D. The lattice 
is broken into 2x2x2 unit cells. Each of these unit cells 
contain eight check operators (one of which is left unde- 
termined) and 33 qubits, nine of which are shared. The 
current distribution over the three walls H,V , and T are 
computed by summing over the bulk configurations con- 
sistent with a given current and excitation configuration. 

There is obviously a computational cost associated to 
summing over the bulk processes of a unit cell. This cost 
grows exponentially with the number of qubits contained 
inside the cell. For this reason, decoding a 2 x 2 x 2 
unit cell involves summing over 26-bit configurations (the 
cell contains 33 qubits and has seven check constraints), 
which is fairly demanding. For this reason, we choose to 
work with smaller unit cells. 

To make the renormalization method for fault- 
tolerance practical, we consider asymmetric decoding. 
The simplest unit cell has dimensions 2x1x1 (see 
Fig. l8|. In this case, the cell contains one magnetic 
fiux operator and renormalizes only one dimension of 
the lattice: £ x € x £ -)► £/2 x £ x ^. For the next 
step, rotate the cell to renormalize another direction, 
e.g. £/2 X £ X € ^ £/2 X €/2 x L Finally, consid- 




FIG. 8: Exploded view of 2 x 1 x 1 unit cell. Qubits live on 
the faces. Qubits 1, 2, 6 and 7 are shared between cells and 
so participate in BP. See Eq. Q for the operator basis. 



ers a second rotation to renormalize the third direction: 
£/2 X £/2 X £ -)> £/2 x (./2 x (./2. For this cell, we choose 
the following operator basis 



So = X1X3X4 
Si = X2X3X5 



To = X3 



Eq = Xi,Xfi 

El = X'sXr 



I/O — X0X3 
Li = X4 (9) 
L2 ~ X5, 



with the same physical interpretation as in the 2D case. 
We have also considered a 2 x 2 x 1 unit cell with the 
following operator basis (see Fig.l9|: 



So = X0X3X5 
Si = X^XaXgXii 

52 = X2X6X8 

53 = X1X4X5 

54 = X7X10X11 

55 = XiXaXr. 



To — X5X9 

ri = X9 

T2 = Xi\ 



Eq = X<)X\2 

El = X11X13 

£-2 = ^5^9^14 

E-j, = X9X15 
i?4 = X11X16 



IV. NUMERICAL RESULTS 

The 2D version of this RG decoding algorithm was 
numerically benchmarked in [T| [5] for Kitaev's toric 
code, and in [16] for the Z^ generalization of the toric 
code. Here, we present numerical results obtained for 
the 3D fault-tolerant case (see also [H]). We consider the 
isotropic case where every qubit is independently subject 
to a bit-flip noise with probability p and likewise measure- 
ments are subject to independent noise that flips their 
outcome with probability p. We use standard Monte 
Carlo techniques to estimate the fault-tolerant storage 
threshold. Our results are shown in Fig.[TO|for the 2x1x1 
cell and and Fig. [TT] for the 2x2x1 cell. Thresholds 
are observed at pth ^ 1.8(2)% and pth ~ 1.9(4)% respec- 
tively: for p < pth, the failure probability of the decoding 
algorithm decreases as the lattice size increases. These 
values should be compared to the 2.9% value obtained 
via PMA [3] with the same error model. 




FIG. 9: Exploded view of a 2 x 2 x 1 unit cell. Qubits 0, 1, 
2, 4, 7, 12, 13, 14, 15 and 16 are shared. See Eq. ([lO| for the 
operator basis. 



0.01 



I/O — X3X9 

I/i ~ ^10 

1/2 ~ XgXii 

(10) 
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FIG. 10; Monte Carlo estimation of the decoding error prob- 
ability as a function of bit-flip channel strength, p using a 
2x1x1 unit cell. A threshold is observed at ~ 1.82% (sam- 
ple size: IC* per point). 



Note that the 2x2x1 unit cell is only compatible 
with lattice sizes that are powers of four. Moreover, due 
to the large size of the unit cell, decoding is relatively 
slow in this case, which limits us to small lattices t — IQ 
and £ — 6A in practice.^ The crossing point of the cor- 
responding two curves gives us little confidence that we 
have correctly identified the threshold. For this reason. 



* The complexity of the RG scheme is proportional to the space- 
time volume of the lattice, while the complexity of PMA scales 
with the cube of this volume. However, the constant factor of 
the RG scheme is exponential with the volume of each unit cell. 
Although this is independent of the lattice size, the constant can 
be quite prohibitive for large unit cells. Note also that RG can be 
straightforwardly parallelized to run in time logarithmic with the 
space-time volume of the lattice, but we have not implemented 
this parallel version. 
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FIG. 11: Monte Carlo estimation of the decoding error prob- 
ability as a function of bit-flip channel strength, p using a 
2x2x1 unit ceU. A threshold is observed at ~ 1.9(4)%. We 
note some finite size effects for ^ = 8, but not for the other 
three curves. Sample size varies from 3 x 10"^ to 10*. 



we also simulated lattice sizes .^ = 8 and £ = 32 using an 
hybrid techniques where the 2x2x1 cell was used until 
the very last step, where a 2 x 1 x 1 call was used. The 
crossing point of all four curves agrees very well. This is 
not surprising since below threshold, we expect the er- 
ror model to flow to a noiseless fixed-point, and therefore 
the failure rate should be largely independent of how de- 
coding is performed at the last few RG iterations — the 
first RG iterations are the critical ones in determining 
the threshold. This observation also gives us confidence 
that RG could handle various lattices shapes by combin- 
ing different unit cell shapes in the appropriate way deep 
in the RG flow. 

One might suspect the threshold to be anisotropic — 
given the asymmetry in the RG, e.g. the direction that 
is renormalized first might exhibit a lower threshold. We 
analyzed the data by looking at the marginal error rate 
in the three different directions and saw no significant 
anisotropy. It is possible that the threshold is insensitive 
to such details, but that they have more subtle effect such 
as leading to different scaling exponents. In both cases, 
better statistics would be needed to give a quantitative 
answer. 



V. CONCLUSION AND OUTLOOK 

We have given a detailed presentation of a renormal- 
ization group algorithm for fault-tolerant decoding of 
topological quantum error-correcting codes supporting 
Abelian anyonic excitations. This extends our previous 
work [3 H] in an essential way, permitting error correc- 
tion in the presence of faulty measurements. We have 
numerically benchmarked this algorithm and found that 
it achieves a fault-tolerant error threshold of nearly 2%, 
in the same ballpark as the other leading techniques. 



A. Relation to other work 

Since the publication of our algorithm [li ^ , there has 
been a number of decoding algorithms proposed for topo- 
logical codes that we now briefly review. 

Sarvepalli and Raussendorf (SR) [H] have conceived 
a RG decoder for topological color codes that resembles 
ours in many ways. To our understanding their algorithm 
is conceptually identical to ours. Their presentation dif- 
fers in one central way. Because some stabilizer genera- 
tors unavoidably overlap with two different cells we were 
forced to share qubits between unit cells, which led to 
inter-cell correlations. Instead of this, SR split those sta- 
bilizer generators into two parts, each supported on a 
unique cell, and assign a binary random variable to the 
value of each part. The sum of these two random vari- 
ables must equal the value of the syndrome associated 
to the stabilizer. These auxiliary binary variables play 
a role analogous to the shared qubits in our description. 
For the color code, their decoder achieves a threshold of 
7.8%, compared to 8.7% achieved by mapping the code 
to multiple copies of KTC and decoding them with our 
RG algorithm [13]. 

Bravyi and Haah (BH) ^20^ have proposed a RG de- 
coder suitable for any topological code supporting lo- 
calized Abelian anyons. It crucially differs from our 
approach by being based on hard decisions, while our 
approach uses soft decisions. In other words, the opti- 
mal recovery is only decided at the very last step of our 
RG iterations. At intermediate iterations, probabilities 
are assigned to various recoveries, but none of the op- 
tions is ever ruled out until the very end. In contrast, 
in the BH scheme, decisions are taken to fuse certain 
pairs of excitations at intermediate iterations of the RG 
scheme. Hard decoders are conceptually simpler, and so 
lend themselves to more rigorous analysis. Indeed, BH 
were able to prove that their decoder achieves a flnite 
threshold, while we can only provide numerical evidences 
for our algorithm. On the other hand, it is well known in 
classical coding theory that soft decoders achieve better 
performances |1T]. In the quantum setting, it has been 
shown that soft decoder can achieve a higher threshold 
and greater noise suppression below threshold |22j . Their 
algorithm achieves a threshold of 6.7%. 

Wootton and Loss (WL) [53] used Monte Carlo sam- 
pling to estimate the sum in Eq. Q, thus directly es- 
timating the probability of each equivalence class of er- 
rors conditioned on the error syndrome. Since Monte 
Carlo is exact within statistical error, given a sufficiently 
large sample, this technique is optimal and consequently 
outperforms all other decoding algorithms. Indeed, they 
achieve a threshold of 18.5%, compared to 16.4% using 
our method with the same noise model. Its main draw- 
back is that it is very slow compared to other methods, 
its runtime scales (morally) exponentially with the lattice 
size. 

Lastly, Fowler, Whiteside and HoUenberg (FWH) [H] 
have implemented a parallelized version of Edmonds' per- 
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feet matching algorithm [TT] (PMA), which was the first 
algorithm used to decode topological code [8j. This im- 
plementation runs in constant average time without any 
performance loss compared to the original PMA. Our un- 
derstanding of this algorithm is that it is of Las Vegas 
type, meaning that its run-time is not pre-determined. 
For instance, in this parallel implementation, it is possi- 
ble that one node of the cluster requires more time than 
other nodes. On a very large lattice, these fluctuations 
could become important, i.e. the probability that at least 
one node takes a time superior or equal to any finite T 
approaches one. Thus, care must be taken in the inter- 
pretation of this constant average runtime. 



B. Extensions 

It is possible to combine these techniques in various 
ways to obtain tradeoffs between runtime and error cor- 
rection. For instance, the RG algorithm of BH can con- 
ceptually be seen as a degradation of our algorithm where 
probabilities on current variables V{1) are rounded up to 
the closest binary distribution 



V'{1) 



1 if I maximizes V{1) 
otherwise 



(11) 



Because of this simplicity, it is much faster than our algo- 
rithm. There exist intermediate degradations that could 
interpolate between these two extreme schemes. For in- 
stance, we could round up the distribution to the closest 
trinary distribution 



V'{1) 



1 if Vil) > 1 
if Vil) < e 
F otherwise 



(12) 



where the flag symbol F is used to signal a potential er- 
ror. Such a scheme was used by Knill [53] in the context 
of concatenated codes, which can be seen as a degrada- 
tion of the scheme of [22] that uses the exact probability 
distribution. One could also consider keeping only the 
first few largest probabilities, and rounding all other to 
zero. 

As we have seen, larger unit cells lead to better er- 
ror correction, but the exact summation Eq. ^ of bulk 
processes inside a unit cell scales exponentially with the 
volume of the cell. One possibility would be to sum the 
bulk processes inside the unit cell only approximately. 
This would enable RG decoding using much larger unit 
cells. For instance, we could use WL's Monte Carlo's 
scheme to estimate this sum. Alternatively, we could 
use tensor-network techniques [26^ to approximate this 
sum. Even without approximations, a transfer matrix 
approach could be used to decrease this complexity from 
exponential in the area of the cell (or volume in 3D) to 
exponential in its linear size (or area in 3D). For the small 
cells we considered here, these more elaborate techniques 
are of no use. 



Lastly, we note that the description of our algorithm 
presented in Sec. |ll] applies equally well to subsystem 
codes [17] that have local stabilizer generators in 2D, such 
as the topological subsystem color codes [T5] (but ex- 
cludes e.g. Bacon-Shor codes [2S]). Indeed, the stabilizer 
generators of these codes reveal excitations that carry 
topological charges and the decoding problem consists 
of inferring the world-line homology of these excitations. 
The main difference is that not all topological charges can 
corrupt the encoded information. Some of the topolog- 
ical charges — that we called gauge charges in |13| — can 
be dragged along a non-trivial cycle without changing the 
ground state of the system. Thus, the current associated 
to these charges does not need to be monitored. Thus, 
Eq. ([6]) should contain an extra sum corresponding these 
harmless processes. We have used this technique for the 
topological subsystem color code in [13 and obtained a 
threshold of 1.95%. 
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Appendix A: Manipulating probabilities over Q" 

In this Appendix we provide some mathematical back- 
ground for manipulating probabilities over the n qubit 
Pauli group Q^\ This should be useful to understand the 
details of Sec. |III| or to implement the RG algorithm. 

Let V{E) be a probability distribution over the n-qubit 
Pauli group t/", e.g. corresponding to a physical noise 
model. Given a generating set {Qi} of t/", we can express 
any E G G"^ as 



E ■ 



4=1 



(Al) 



where Xi £ {0, 1}. This allows us to interpret ViE) as 
a distribution over 2n binary variables V{xi, ...,X2n) = 
ViE — Y[i=i QT)- Standard Bayesian calculus can then 
be used to define marginal distributions, conditional dis- 
tributions, etc. For instance, the marginal distribu- 
tion over xi, X2, and x^ is given by V{xi,X2,x^) = 
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Sx4,...i:2„ ^(^i'^2,- • •2:2n)- The probability of Xi and 
X2 conditioned on x^ is given by V(xi,X2\x-i) — 
V{xi,X2,X3)/V{x3). These probabilities implicitly de- 
pend on a basis choice {Qt}, and we can perform such 
manipulations for any basis of Q"". 

These definitions extend straightforwardly to more 
variables. With the isomorphism Eq. (Al), we can re- 
label these probabilities ^^(Qi, Q2, Q3) = 'P{xi,X2,X3), 
and so forth. 

We can create coarse grained variables associated to 
subgroups of the Pauli group. For instance, let K, — 
{Qi,Q2,Q3} and T = (QiiQb) be two subgroups of 
5". An element K oi K, can be decomposed a,s K = 
Qi^ QT Qs^ 1 and similarly an element T of T can be de- 
composed as -ftT = Q^^Q^^. The joint, marginal, and 
conditional probabilities can then be defined in a natural 
way 

ViK, T) - V{KT) = V{xi,X2,X3,Xi, x^) (A2) 

r{K)=r{xuX2,X3) (A3) 

r{T) = r{xi,x5) (A4) 

r{K\T)^Vixi,X2,X3\xi,X5). (A5) 



These are the formal definitions behind Eqs. ( 4|6|7 ). 

Lastly, we can convert any of these probabilities — ^joint, 
marginal, and conditional — to a different basis. For in- 
stance, let (Qi,Q2,Q3) be a different generating set for 



JC. We can express these generators in terms of the pre- 
vious ones Q[ = Ylj=i 2 3 Q^'^ with yij G {0, 1}. Suppose 
that we have computed V{K\T) = V{xi,X2,X3\x4,X5) 
using the basis {Qi}, and now wish to compute 'P{K\T) 



2"-Q'3'. Since 




K= UiU Qrr 


(A6) 


i=l,2,3 j = l,2,3 




= IT Qf--'^"^% 


(A7) 



J = l,2,3 

we see that 'P{K\T) = 'P{zi,Z2,Z3\x4,x^) — 
'P{xi,X2,X3\xi,X5) for Xj = J2i=i,2,3y'-3^i- Thcsc prob- 
abilities can then be used to compute marginals over a 
subgroup of /C specified in terms of the primed genera- 
tors. For instance, for F G {Q'i,Q2) we have r{F\T) = 
V{zi, Z2\x4,X5). Thus, we see the usefulness of perform- 
ing basis changes: it is used to adapt the probability to 
the particular subgroup we are interested in. 

We will be using this type of manipulation in the spe- 
cial case where the basis {Q'A actually corresponds to 
the basis of single qubit Pauli operators {Xi, Zi}. In that 
case, for ii^ = Hi ■^i'' ^i ^^ ^^^^ be using the special no- 
tation K\q to represent Xq'' Zq , i.e. the Pauli operator 
on qubit q in K. These are the formal definitions behind 



many mathematical expressions of Subsection III D 
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